#include "TH2.h"
#include <iostream>
#include <cmath>

float energyDensity(float,float,float,float,float,float);

TH2D* coneDepot(float energy,float polarAngle, float azimuthalAngle)
{
	TH2D *h2Cone = new TH2D("h2Cone","Cumulative Energy Deposits;X (cm);Y (cm)",40,-50,50,40,-50,50);

	float d = 0.001;

	for(int i = 27;i<48;i++)
	{
		for(int j = 18;j<29;j++)
		{
			float iCenter = (i-24+.5)*0.025;
			float jCenter = (j-24+.5)*0.025;

			for(float iFine=(-0.025+d)/2;iFine<.025/2;iFine=iFine+d)
			{
				for(float jFine=(-0.025+d)/2;jFine<.025/2;jFine=jFine+d)
				{
					for(float kFine=-d/2;kFine>-.5;kFine=kFine-d)
					{
						float iAbs = iCenter+iFine;
						float jAbs = jCenter+jFine;
						float kAbs = kFine;

						float deposit = energyDensity(energy,polarAngle,azimuthalAngle,iAbs,jAbs,kAbs)*pow(d,3);
						if(deposit>0)
						{
							h2Cone->Fill(iCenter*100,jCenter*100,deposit);
						}
					}
				}
			}
		}
	}

	return h2Cone;
}

float energyDensity(float energy,float polarAngle,float azimuthalAngle,float i,float j,float k)
{
	float C = 0.0186/log(2);
	float energyC = 0.07;
	float depth = C*log(energy/energyC);
	float radius = 0.038;

	float distanceFromStart = pow((pow((i-(12+.5)*0.025),2)+pow(j-(0+.5)*0.025,2)+pow(k,2)),.5);
	float newAzimuthal = acos(((i-(12+.5)*0.025)*sin(polarAngle)*cos(azimuthalAngle)+(j-(0+.5)*0.025)*sin(polarAngle)*sin(azimuthalAngle)+k*cos(polarAngle))/distanceFromStart);
	float maxAzimuthal = atan(radius/depth);

	if(newAzimuthal < maxAzimuthal)
	{
		if(distanceFromStart < depth/cos(newAzimuthal))
		{
			float C2 = 12*energy/(pow(depth*radius,2)*3.14159265);
			return (depth-cos(newAzimuthal)*distanceFromStart)*C2;
		}
	}

	return 0.0;
}
